GHOST FORCE INFLUENCE OF A QUASICONTINUUM 
METHOD IN TWO DIMENSION 



JINGRUN CHEN AND PINGBING MING 

Abstract. We derive an analytical expression for the solution of a two-dimensional 
quasicontinuum method with a planar interface. The expression is used to 
prove that the ghost force may lead to a finite size error for the gradient of the 
solution. We estimate the width of the intcrfacial layer induced by the ghost 
force is of 0(y/e) with e the equilibrium bond length, which is much wider 
than that of the one-dimensional problem. 



1. Introduction 

Multiscale methods have been developed to simulate mechanical behaviors of 
solids for several decades [14]. Combination of models at different scales greatly 
enhances the dimension of problems that computers can deal with. However, prob- 
lems regarding the consistency, stability and convergence of the multiscale methods 
may arise from the coupling [3J. Taking the quasicontinuum (QC) method [171 110] 
for example, one of the main issues is the so called ghost force problem |16j , which 
is the artificial non-zero force that the atoms experience at the equilibrium state. 
In the language of numerical analysis, the scheme lacks consistency at the interface 
between the atomistic region and the continuum region [4j . For the one-dimensional 
problem, it has been shown in |15j that the ghost force may lead to a finite size 
error for the gradient of the solution. This generates an interfacial layer with width 
C(e|ln£|), out of which the error for the gradient of the solution is of 0(e). 

To understand the influence of the ghost force for high dimensional problems, we 
study a two-dimensional triangular lattice model with a QC approximation. This 
QC method couples the Cauchy-Born elasticity model and the atomistic model with 
a planar interface. Numerical results show that the ghost force may lead to a finite 
size error for the gradient of the solution as in the one-dimensional problem. The 
error profile exhibits a layer-like structure. It seems that the width of the intcrfacial 
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layer induced by the ghost force is of 0(y/e), which is much wider than that of the 
one-dimensional problem. 

To further characterize the influence of the ghost force, we introduce a square 
lattice model with a QC approximation. Compared to the triangular lattice model, 
this model can be solved analytically and the error profile exhibits a clear layer-like 
structure. Based on the analytical solution, we prove the error committed by the 
ghost force for the gradient of the solution is 0(1) and the width of the induced 
interfacial layer is of 0(y/e), which are also confirmed by the numerical results. 

The paper is organized as follows. Numerical results for the triangular model 
and the square lattice model with QC approximations are presented in §[2] and §[3j 
respectively. We derive an analytical expression of the solution of the square lattice 
model with a QC approximation in § 0J The main results of the paper are proved 
in§[5l 

2. A QUASICONTINUUM METHOD FOR TRIANGULAR LATTICE 

2.1. Atomistic and continuum models. We consider the triangular lattice L, 
which can be written as 

L = { a; G K 2 | x = ma\ + na2, m,n £ Z} 

with the basis vectors a\ = (1, 0), a 2 = (1/2, y/3/2). Define the unit cell of L as 

r = {ieM 2 I x = aai + c 2 a 2 , -1/2 < ci,c 2 < 1/2}. 

We shall consider lattice system eh inside the domain f2, and fl £ = ft n eh, where 
e is the equilibrium bond length. Assume that the atoms are interacted with the 
potential function V, which is usually a highly nonlinear function, e.g., the Lennard- 
Jones potential [12]. Denote by S\ and S 2 the first and the second neighborhood 
interaction ranges; see Figure [T] In particular, we have 

51 = U i=1 Sj = {di, 02, —a\ + a 2l —ai,—a2) 01 — a 2 }, 

5 2 = L>]l 7 Si = {ai + a 2 , -ai + 2a 2 , -2a% + a 2 , -a\ - a 2 ,ai - 2a 2 ,2ai - a 2 }. 

For (i e Z 2 , the translation operator is defined for any lattice function z : 

L — s- R 2 as 

(T^z)(x) = z(x + e/iiai + E[i 2 a 2 ) for x € L. 
We define the forward and backward discrete gradient operators as 

D+ = £- 1 (T»-I) and D~ = e~ 1 (J - T £ " ) , 

where s = [iiai + fi 2 a 2 and / is the identity operator. We shall also use the short- 
hand Dz = (Dfz,D£z) = (D+z,D+z). 
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FIGURE 1 . The first and second neighborhood interaction ranges of the 
triangular lattice; <Si = {si, ■ ■ ■ , S&} and S2 = {S7, ■ ■ ■ , S12}. 

Consider an atomic system posed on f2 e . The total energy is given by 

(2.1) £* a f = ^E E V(\D+z(x)\), 

where V is a potential function. In this paper, we only consider the pairwise 
potential function, and leave the discussion on the more general potential functions 
in future publication. The atomistic problem is to minimize the total energy subject 
to certain boundary conditions that will be specified later on. 

Next we turn to the Cauchy-Born elasticity model [TJ [5J [5J [7] • Given a 2 by 2 
matrix A, the stored energy density function is given by 

W cb (A) = -L E V(\s-A\), 

sG5iU5 2 

where i?o is the area of the unit cell and t?o = v3e 2 /2. The stored energy function 
is defined by 

= [ W ch (Vz(x))dx. 
Jn 

The continuum problem is to minimize the stored energy function subject to certain 
boundary conditions. We employ the standard P\ Lagrange finite element to ap- 
proximate the Cauchy-Born elasticity model with the lattice L as the triangulation. 
The approximate stored energy function is 

(2.2) EX £ = i E E i\ D tA*)\) + V (UK + Dt i+ Mx)\)) . 

a;£fi e i— 1 

One can see E^ e reproduces the atomistic energy ; cf. (|2.1[) , if only the nearest 
neighborhood interaction is considered. 

We study the quasicontinuum method [17]. Let e — 1/(2M). We assume that 
the interface between the continuum model and the atomistic model is x\ = as 
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shown in Figure [2j The total energy of the QC method is 

E T = l E Y,( v (\ D i z w\)+ v {\( D t+ D t i+ >(*)\)) 

" iri<-2e i=l 

+ \H {E(^(i^^(^)i) + ^(i(^+^+ +1 M a; )i))+ y (i^ + 1 ^^)i) 

xi = — e i=l 

^ E { E ^ i\ D t^)\) + E (V (\ D tA*)\) + v {\( D t + D f i+ M*)\ 

— n £z,c, 4—0 



2 

xi=0 v seSi i=2 

+ + V(\Df iiZ (x)\) + V(\Df i2 z(x)\) 

^E{ E V(\D+z(x)\)-±V(\D+z(z)\)\ 

J E E ^(i^)D- 



Ki>2e sS5iU5 2 

The force at atom x is defined by 

Since we only concern the influence of the ghost force, following |15j . we assume 
that the potential function is a harmonic function, i.e., 

V(r) = \r\ 

Without taking into account the external force, we write the equilibrium equations 
for the QC approximation as 

F qc [z](x) = 

with 

12 

J"q C [^](x) = — 12z(x) + z(x + esj), x e 0. s ,xi < — 2e, 

»=1 
6 

J" qc [z](x) = — 24z(x) + 4^ z(x + esj), x e xi > 2e. 

i=i 

For a; = (— e,x 2 ), 

49 1 

^qeWw = — ^z{x) +42^z(i + £Sj) + -z(x + esi 2 ). 

i=l 

For x = (0, x 2 ), 



2 l 

+ 4[,z(x + es 3 ) + z(x + es 4 )] + [z(x + es 7 ) + z(x + esu) + z(x + esia)]. 



•Fq C [z](:c) = -18z(x) + [z(x + esi) + z(x + es 6 )] + -[z(x + es 2 ) + z(x + es 5 )} 
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m=-M m =° m=M 

Continuum region Atomistic region 

FIGURE 2. Schematic picture of fi e . 



For x = (e, x^) 



23 ^ 

J" qc [z](a;) = -—z{ps) + ^z(x + esi) ~ -z(x + es g ). 

i=l 



At the equilibrium state, we evaluate J-" qc at z(x) 



T qc [x](x) 



( (-3e/4,V3e/4), 
(3e/2,-VSs/2), 
(-3e/4, V3e/4), 
(0,0), 



2 

x to get 

Xl = -£, 

xi =0, 



if 
if 

if Xl = £ 
otherwise 



The above equations imply Zx(x) = — y/3z2(x). Therefore, we only study the first 
component of z(x) and neglect the subscript if no confusion will occur. The error of 
the displacement induced by the ghost force, denoted by y(x) = z(x) — x, satisfies 



(2.3) 
with 



Fqc[y](x) = T qc [z](x) - T qc [x](x) = -J" qc [x](a;) = f(x) 



3e/4, 
-3e/2, 
3e/4, 




if 
if 
if 



Xl 
Xl 
Xl 



0. 



otherwise. 



Boundary conditions need to be supplemented to close the system of equilibrium 
equations. Two types of boundary conditions will be considered in this paper. One 
is the periodic boundary condition in X2 direction while homogeneous Dirichlct 
boundary condition in x\ direction, which will be called periodic boundary condi- 
tion if no confusion will occur. The other is the homogeneous Dirichlct boundary 
conditions in both x\ and X2 directions. 

In what follows, we shall use a conventional notation y(m,n) = y(x) with x = 
e{ma\ + 7102). 
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2.2. Periodic boundary condition. The periodic boundary condition can be 
written as 

y(-Af,n) = 0, n = -N,---,N, 

y(m,n) = y(m,2N + n), m = —M, ■■■ ,M,n = -N, ■■■ ,N, 

y (M, n) = y(M + 1, n) = 0, n = -N, ■■■ , N. 

Observe that the solution of ()2.3|) with the periodic boundary condition takes a 
special form: 

y(x) = cj/(xi), 

namely, the solution is constant along xi direction. Based on this observation, we 
conclude that the QC approximation with the periodic boundary condition reduces 
to a one-dimensional problem. The equilibrium equations satisfied by yix\) arc as 
follows. In the continuum region, i.e., m = — M + 1, • ■ ■ , —2, 

8y(m + 1) - 16y(m) + 8y(m - 1) = 0, 

and in the atomistic region, i.e., m = 2, • ■ ■ , M — 1, 

y(rn + 2) + 4y(m + 1) - 10y(m) + Ay(m - 1) + y(m - 2) = 0. 

The equations in the interfacial region are 

/I 33 3 

-y(l) + 8y(0)- Y y(-l) + 8y(-2) = -S, 

! /(2)+4y(l)-13y(0) + 8y(-l) = -|e, 

19 13 

^ 2/(3) + 4y(2) - — y(l) + 4y(0) + -y(-l) = -£. 

The boundary condition is 

y(-M) = y(M) = y(M + 1) = 0. 

The above equilibrium equations can also be obtained by considering a one- 
dimensional chain interacted with the following harmonic potential: 

(2.4) v({ y }) = y \y*-y 3 \ 2 + Y E \»-m\ 2 

|i-j'l=2 

with ki = 4 and /c 2 = 1. This is the model studied in [2]. 

2.3. Dirichlet boundary condition. The Dirichlct boundary condition can be 
written as 

y(-M,n) = y(M,n) =y(M + l,n) = 0, n = -N,---,N, 

y(m,-N) = y(m,N) = 0, m = —M, ■■ ■ , — 1, 

j/(?n, -AT - 1) = j/(m, —AT) = y ( m , iV) = ?/(m, iV + 1) = 0, m = 0, ■ ■ • , M. 

This gives an essentially two-dimensional model. We show profiles of discrete gra- 
dients Df t y and Df 2 y in Figure [3] with M = N = 20. The profile of D^y is similar 
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to that of the one-dimensional problem. However, if one zooms in the interface, 
there are some differences. To highlight the interface, for i = 1, 2, we define 

Ei = {xen £ | \D Si y{x)\ >c e}, 

and let \Ei be their characteristic functions. We plot \Ei (%) in Figure 0] with 



M = N = 640 and Cq = 0.04. In Figure 4(a) the width of the interface near 



the boundary is much wider than that in the interior domain. No interfacial layer 
for D+ y is observed in Figure |4(b)| Therefore, we only measure the width of the 
interfacial layer for D^y, which is defined by 

max < |ai|erargmax m {x = ma\E + na^e £ Q £ \xi(x) = 1} 

— argmin TO {x = ma\£ + na 2 e G tt £ \xi(x) = I}]}- 

Numerical results in Table [1] imply that the width of the interfacial layer scales 
0(y/e), while the one-dimensional problem has an interfacial layer with C(e|lne|) 
width PUJH]. 



0.2 
0.1 

-0.1 



-0.2-J 
1 





(a) D+y 



(b) D+y 



FIGURE 3. Profiles of the discrete gradients for triangular lattice with 
M = N — 20 under Dirichlet boundary condition. 



e (5 x 10" 2 ) 2° 2- 1 2- 2 2- 3 2" 4 Rate 
Layer width(I0- 1 ) 3.5 2^3 L5 0.94 0.63 0.52 0.57 

TABLE 1. Width of the interfacial layer versus the equilibrium bond 
length e for the triangular lattice with Dirichlet boundary condition. 



One may doubt the above result could be caused by the boundary condition 
instead of the ghost force. To clarify this issue, we enlarge the continuum region to 
weaken the influence of the boundary condition. We set N = 3M and use different 
equilibrium equations for different regions as in Figure [SJ Roughly speaking, along 
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I , , 1 I — 1 I 1 1 1 L 

0.5 1 0.5 1 

(a) XEi(x) (b) XE 2 (x) 

FIGURE 4. Profiles of the characteristic functions for triangular lattice. 



K1 Continuum region 
m=N „ : 




Continuum region 



FIGURE 5. Padding technique to remove the boundary effect. 

xi direction, the original QC method is employed in the interior region, and the 
continuum model is padded in the outer region. 

We plot D^y, -D+y, and their characteristic functions in Figure [5] and Figure [7] 
They are similar to those in Figure [3] and Figure [4] 

Table [2] shows the width of the interfacial layer in terms of e. Numerical results 
suggest that the width of the interfacial layer is still of 0{yfe), which is consistent 
with that in Table [T] Therefore, we conclude that the interface is caused by the 
ghost force instead of the boundary condition. 



£ (itr 2 ) 


1.7 0.83 0.41 


0.21 


0.10 


Rate 


Layer width(10~ 2 ) 


9.2 8.3 4.8 


3.3 


2.3 


0.53 



TABLE 2. Width of the interfacial layer versus the equilibrium bond 
length e by removing the effect of the boundary condition. 

The QC approximation discussed in this section is quite realistic except the 
potential function as shown in j!5) . which however seems enough to characterize 
the influence of the ghost force. Unfortunately, it does not seem easy to solve this 
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(a) D+y 



(b) D+y 



FIGURE 6. Profiles of the discrete gradients for triangular lattice with 
M = 20 and padding technique under Dirichlet boundary condition. 




(a) X Bl O) (b) XE 2 (x) 

FIGURE 7. Profiles of the characteristic functions for triangular lattice 
with padding technique. 

model analytically as we have done in |15j for the one-dimensional problem. In 
next section, we introduce a new QC method that can be solved analytically. We 
shall prove that this new QC method does capture the main feature of the ghost 
force for the triangular lattice model with a planar interface. 

3. Square lattice model 
We consider the square lattice with the harmonic potential. Compared to the 
standard interaction range of the square lattice, we assume a special interaction 
range as shown in Figure [8] (Left). Namely, the first and second neighborhood 
interactions in x\ direction, and the first neighborhood interaction in x-i direction 
are taken into account. This seemingly strange selection may be obtained from a 
rotated triangular lattice as in Figured] (Right). If we condense the interaction of 
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FIGURE 8. Interaction ranges. Left: First and second neighborhood 
interactions in xi direction and first neighborhood interaction in X2 di- 
rection for square lattice; Right: First and second neighborhood inter- 
actions for triangular lattice in a rotated coordinate. 

the atoms 9, 4, 5, 11 into one atom, and the atoms 8, 2, 1, 12 into another, then we 
obtain a square lattice model with the special interaction range described above, 
which may be the underlying reason why it can be regarded as a surrogate model. 
We shall show in the next two sections that this model not only captures the main 
features of the triangular lattice model with the QC approximation as shown in 
Figure [21 but also lends itself theoretically tractable. 

Proceeding along the same line that leads to (|2.3p , wc obtain the equilibrium 
equations for the error y(a;)Q In the continuum region, i.e., m = —M, • • • , —2 and 
n = —N, ■ ■ ■ ,N, 

12yi(m, n) — yi(m, n — 1) — yi{m, n + 1) — 5yi(m — 1, n) — 5y,:(m + 1, n) 

(3.1) =0, i = 1,2, 

and in the atomistic region, i.e., m = 2, • • • , M and n = —N, ■ ■ ■ ,N, 

6^(771,™) - y t (m,n- 1) - y i (m 1 n + 1) - y^m - l,n) - y t (m + l,n) 

(3.2) -yi(m-2,ri)-yi(m + 2,ri) = 0, i = l,2. 

The interface between the continuum model and the atomistic model is the line 
m = as shown in Figure [9j and M is assumed to be even for simplicity. The 
equilibrium equations for the layers m = — 1, and 1 arc as follows. 

For layer m = — 1 and n = —N, ■ ■ ■ ,N, 
25 

—yi(-l,n) - yi(-l,n- 1) - yi(-l,n+ 1) - 5y.;(-2,n) - 5^(0, n) 

(3-3) -\vi(l,n) = fi, i = l,2, 

where fi = —e and ji = 0. 



^We actually multiply —1 on both sides of 12.31 1. 
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n=N 



m=-M m=0 

Continuum region Atomistic region 



FIGURE 9. Square lattice with a planar interface. 



For layer m = and n = —N, ■ ■ ■ ,N, 

9y t (0,n) - j/j(0,n- 1) - y l {0,n + 1) - 5y l {-l,n) - y,-(l,n) 
(3-4) -Wi(2, ») = /*, i = 1,2, 

where /i = 2e and /2 = 0. 

For layer m = 1 and n = —TV, • • • , TV, 

y2/i(l,n) -j/j(l,n- 1) -2/i(l,n + l)-j/i(0,n)-yj(2,n) 

(3.5) - -yi(-l,n) -2/j(3,n) = /j, i = 1, 2, 

where f\ = — e and ji = 0. 
Observe that 

2/2 = 0. 

Therefore, we only consider y\ and omit the subscript from now on. 

Wc first impose the Dirichlet boundary condition in the x\ direction, and the 
periodic boundary condition in the X2 direction as 



(3.6) 



y{-M, n) = y(M, n) = y(M + 1, n) = 0, n = -TV, • • • , N, 
y(m, n) = y(m, n + 2N), m = —M, ■ ■■ ,M,n = -N, ■ ■■ , N. 



Similar to the triangular lattice model, it is easy to check that the square lattice 
model reduces to a onc-dimcnsional chain model with the following equilibrium 
equations and boundary condition. 



5y(m + 1) - 10y(m) + 5y(m - 1) = 0, m= -M + 1, • • • , -2, 
y(m + 2) + y(m + 1) — Ay{m) + y(m — 1) + y(m — 2) = 0, m = 2, • • ■ , M — 1. 
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The equations for the interface are 
1 



21 



2/(1) + 5y(0) - yZ/M) + 5y(-2) = e, 
y(2)+y(l)-7y(0) + 5y(-l) = -2e, 



j/(3) + 2/(2) - |z/(l) + 2/(0) + iy(-l) = e. 

The boundary condition is 

y(-M) = y(M) = y(M + 1) = 0. 

It is clear that the above model is exactly the same as that has been studied 
in [15], which can also be obtained from the one-dimensional model (|2.4p with 
ki = k 2 = 1. 

Therefore, we impose the Dirichlet boundary condition in both x\ and x 2 direc- 
tions as follows. 



(3.7) 



2/(-M, n) = y{M, n) = y(M + 1, n) = 0, n = -N, ■ ■ • , iV, 
2/(m, -iV) = y(m, N) =0, m = —M, ■■■ ,M. 



Choosing M = N = 20, we show profiles of the discrete gradients Df y and 
Df 2 y in Figure [TU1 The feature is similar to that of the triangular lattice model. To 
highlight the interface, we plot the characteristic functions XE t and \e 2 m Figurc fTTI 
with M = N = 320 and c = 0.04. 




We report the width of the intcrfacial layer in Table |3l It is clear that the width 
of the interfacial layer is of 0(y/e), which is consistent with that of the triangular 
lattice. In § [5j we shall prove this fact. 
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0.6 0.8 



(a) XEiO) 



(b) XE 2 (x) 



FIGURE 11. Profiles of the characteristic functions for square lattice. 



e (5 x 1CT 2 ) 



2° 2- 



1-2 o-3 o-4 



Rate 



Layer width (1CT 1 ) 8 5.5 3.63 2.56 1.84 0.53 

TABLE 3. Width of interfacial layer versus the equilibrium bond length 
e for square lattice. 



4. Exact solution for the square lattice model 



To find the exact solution of the QC approximation (|3.1[) - (|3.5|) with Dirichlet 
boundary condition (|3. 7[) , we follow the approach in [T5] : firstly, we find the general 
expression for the solution of the continuum equation and the atomistic equation by 
separation of variables ansatz, with certain unspecified constants; secondly, we use 
the equations around the interface to determine these constants. The next lemma 
gives the general expression of the solution. 

Lemma 4.1. For m = —M, ■ ■ ■ , — 1 and n = —N, • • ■ , N , we have 

2N-1 



(4.1) 

where 



(m,ri) = ^2 a k smh[(M + m)a k ] sin -^-(N + ti), 



y 



k=l 



2N 



COsh Qffe = 1 



Aj. = 2 sin^ 



kit 
4JV' 



For m = 0, • • ■ , M and n = —N, • ■ • , N, we have 



2W-1 



(4.2) y(m, 
where 



)= E (^(- 1 ) m ^(7 fc ,4) + c fe / m .(7 fc ,4))sin|^(7V + n), 



fc=i 



, 1 + \/25 + 8A fc -1 + V25 + 8A fe 
(4.3) cosh7 fc = , cosh 4 = 
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(4.4) 



sinh[(Af + 1 — 771)7] + sinh[(M — 771)7] cosh 6 



- -1 



cosh 7 + cosh 6 
iC0sh[(7l/ — m)S] sinli7 
cosh 7 + cosh 5 



. ^(7^)=^^, 7)- 
The coefficients bk and Ck are parameters to be determined; See; cf, (|4.13[) . 



Proof. By separation of variables, we get (|4. 1|) . 

The explicit expression for the solution of the atomistic model can also be 
obtained by separation of variables ansatz. Substituting y(m,n) = f(m)g(n) 
into (|3.2[) . we get 

2 1 
E I/( m + - /(*(«) + E + *') " 5(»)]/(m) = 0. 

i=-2 i=— 1 

By (|3.7p . we have 

ff(-iV) = ff(iV) = and f(M) = f(M + 1) = 0. 
We write the above equation as 

SL- 2 [/(^ + 0-/M] , El-Mn + l )-9(n)} =Q 

For A £ R, we get 

J g(n + l) + (2X-2)g(n)+g{n-l)=0, 

\ f(m + 2) + f(m + 1) - (4 + 2A)/(m) + /(m - 1) + /(m - 2) = 0. 

Using the boundary condition for g, i.e., g(N) = g(—N) = 0, we have, for any 

c e M, 

g(n) = csin — (n + N) and A = 2 sin 2 — . 
The characteristic equation for /(m) is: 

i 2 + ?T 2 + f + tT 1 - 2(A + 2) = 0. 
Denote the roots of the above equation by t±, ■ ■ ■ , £4. It is clear that 
h 

with 



-e>, t 2 = -e-\ t a = e s , i 4 = e~ 5 , 



2 cosh7 = — si, Si = 



-1 - V25 + 8A 



2 cosh 5 = S2 , S2 = 



-1 + V25 + 8A 
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This leads to 

/(to) = a(-l) m sinh[(M - 771)7] + b(-l) m cosh[(M - 771)7] 
+ csinh[(M - m)S] + dcosh[(M - m)S] 
with constants a, b, c and d that will be determined by the following conditions 

f(M) = f(M + 1) = 0. 
Since M is even, by f(M) = 0, we obtain 

b = -d. 

By f(M + 1) = 0, we obtain 

a sinh 7 — c sinh 6 — 6(cosh 7 + cos S) = 0. 

Therefore, 

a sinh 7 — c sinh S 
cosh 7 + cosh 5 

We write /(to) as 

/(to) = a(— l) m sinh[(M — 771)7] + csinh[(M — to)<5] 

a sinh 7 — c sinh (5 r , N m , . , . , , . . . „. 1 

+ ^ — -i (-1) cosh (M - 771)7 - cosh M - m)S] [. 

cosh7 + cosh(5 L v ' LV ' ' LV ; J i 

It is easy to rewrite /(to) into a symmetrical form 

f(m) = a(-l) m F m ( 7 ,6) + cF m (5, 1 ), 
where -F m (7, 5) is given in (|4.4[) . this gives (|4.3[) . □ 

Remark 4.2. XTie exact solution based on the series expansion is common in finite 
difference. We refer to |13| and [5] for a thorough discussion. 

Next we use the interfacial equations p.3[) - (|3.5p to determine the coefficients 
, 6fe and . Denote 

Vcb(jn, Ti) = y(m, n) — M < m < — 1 , — N < n < N, 
Vat(jn, n) = y(m, n) < m < M, —N < n < N. 

Though ?/ c b and 7/ a t are only defined piecewisely, they actually admit a trivial 
extension for any (to, n) €E 1? as 

(4.5) J r e [y ch ]{m,n) = 0, _F at [y a t] (to, n) = 0, for (m,n)eZ 2 , 

where 

T e [y] (to, n ) = 12y(m, 77) — t/(to, 77 — 1) — t/(to, n + 1) — 5y(m — l,n) — 5y(m + 1,7),), 
■FatM^n) = 67/(777,77) - y(m,n - 1) - 7/(777., 77 + 1) - 7/(777 - l,ri) - y(m + l,ri) 

— 7/(777. — 2, 77) — 7/(777. + 2, n). 
This observation is crucial to simplify the equations around the interface. 
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The equation for m = — 1 changes to 

F e [Vcb{i,n)] + 5(?/cb - 2/at)(0,n) + ^(y c h{l,n) - y a t(l,n) 

Using 

Fe[Vcb(l,n)] = 0, 

we have 
(4.6) 

Proceeding in the same fashion, we get 
(4.7) 

(4.8) 3j/ at (0, n) + y at (l, n) + j/ at (2, n) - 5y cb (l, n) = 2e 



-e. 



5(y c b - y at )(0,n) + -(2/ cb (l 



In what follows, we use the above simplified intcrfacial equations and the repre- 
sentation formulas to determine a k , bk and c k - 
Subtracting (|4.7j) from (|4.6|) . we obtain 

%cb(0, n) +y ch {l,n) = 5y at (0,n) +y at (l,n). 

Substituting (|4. 1[) and (|4.2[) into the above equation, we get 

& fe (5F - *i)(7fc,4) + c fe (5/o + /i)(7fe,5 fe ) 



(4.9) 



5sinh[Mafc] + sinh[(A/ - l)a fe ] 
Substituting (|4TTj) and P~2"|) into (ET71) . we get 

2JV-1 



4 sin 



2N 



(n + N) 



= -2e 



with 



4 = - sinh[(M - l)afc]a fc + (-2*1 + Fi)(7 fe , + (2/j - /i)( 7fe , S k )c k . 
Using the discrete Fourier transform, we get 



. 2 x (-2s) 2 ^ 1 . fcTTj 

4 = > sm 

27V - 1 + 1 ^ 2N 



-COt— , lf fclS0dd, 

0, if fc is even. 



This leads to 

(4.10) 

where 

(4.11) 



Pkbk +PkCk = 4, 

f P fe = [-2F T + F 1 - p k (-F t + 5F )]( 7fc , 5 k ), 

Pk = [2/i - fx - Pk(h + 5/ )](7fc, 5 k ), 
sinh[(M - l)ajfe] 



5sinh[Ma fe ] + sinh[(M - l)a k ] ' 
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Using (|4.7j) and (|4.8j) to eliminate y c b(l,n), we obtain 

3?/at(0,n) - 9y at (l,n) + y at (2,n) + 5y at (l,n) = 12e. 
The coefficients bk and Ck satisfy 
(4.12) R k b k +r k c k = -U k , 

where 

J J? fc = (3F + 9Ff +F- 2 - 5Fi)(rfk,S k ), 
\ r fe = (3/ -9/i + /2 + 5/i)( 7fcl 4). 

To solve the linear system (|4. 10[) and (|4. 12[) . we need to check whether P k rk — PkRk 
is nonzero for all k. We shall prove in Lemma 15.61 that this is indeed the case. 
Therefore, we may solve (|4.10[) and (|4.12j) to obtain 

r k + 6 Pk n Rk + 6Pk 



(4.13) 



bk 



Ck = -- 



Pk rk ~ Pk Rk Pk Tk-PkRk 

Substituting the above equation into (|4.9I) . we get 



Ik. 



= (r k + 6pk)(5F - Fi)( 7fc , S k ) - (R k + 6P fc )(5/ + fi)(-y k ,6 k ) 
ak (5 S mh[Ma k }+smh[(M-l)ak})(P k rk-PkRk) *" 

To sum up, we obtain the representation formulas for the solution of the QC 
approximation by specifying the parameters a k , b k and Ck in Lemma 14.11 



Theorem 4.3. Let y be the solution of (j3T|) - (|3l)j) . Then for m = -M, • 

Q k sinh[(Af + m)a k ] 



1 



and n = —N, ■ ■ ■ , N , 



2N-1 



y 



(m,n) = -- J2 



N ^ P k r k - R k Pk sinh[(M - l)a k ] 

k odd 



Pk 



(4.14) 



kn 

x cot sm 

AN 



(n + N) 

2N y ' 



where p k is given in (|4. 1 1)) . 

For m = 0, • • • , A/ and n = -N, ■■■ ,N, 



(4.15) 



y(m,n) 



2s 
"iV 



2JV-1 

E 

*;=i 

odd 



Qm.fc , kn . 

cot sm 

P k r k - RkPk AN 



klT 

2N 



(n + N) 



where 

P^k - RkPk = 6(8p fc - 1) 
+ (2~Pk) 



F -Ff 
fo fi 

-Fi F- 2 
fl fa 



+ (25 Pfe - 3) 
+ ($Pk - 1) 



-F 1 


Fo 






F-2 


h 


fo 




fx 


fa 


-F 1 

h 


-Fi 
h 




5p k 


F - 
fo 
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12 


F 


-Fi 


+ 5 


F 


F-2 


+ F 1 




fo 


h 




fo 


h 


-h 



-Fi F- 2 



h 



Fx 

h - h 



Q m ,k = 3(1 - Wpk) 
+ 3(l-2 Pfe ) 



■l) m F m 

fill 

-l) m F„ 

frn 



F 
fo 



-Fi 
fi 



{-l) m F„ 

.fill 



F-2 

h 



fi 



As an immediate consequence of the above theorem, the solution is symmetrical 
with respect to n = 0, i.e., 

(4.16) y(m,n) ~ y(m, —n), 

which can be easily verified from the representation formulas (|4.14j) and (|4.15[) . 
5. Estimate of the QC Solution for the square lattice model 
The main result of this section is the following pointwise estimate of the solution. 



Theorem 5.1. Let y be the solution of (|3.1j) - (|3.5[) . Then 



(5.1) \Dy(m,n)\ < —r exp 



(5.2) \Dy(m,n)\<C 




exp 



2m 



157V 



, m > 0. 



A direct consequence of the above theorem is the estimate of the width of the 
interfacial layer, that is, the region beyond which \Dy\ = 0(e). 

Corollary 5.2. The width of the interfacial layer is 0(\fe). 

We exploit the explicit expression of the solution in Theorem 14.31 to prove The- 
orem l5.il Note that the terms Puru — RkPk, Qk and Q m ,k consist of the terms like 
(—l) m F m f n — (—l) n F n f m for different integers m and n. The asymptotical behav- 
ior of such terms will be given in a series of lemmas, i.e., Lemma 15.61 Lemma 15.71 
and Lemma [5.8l We begin with certain elementary estimates that will be frequently 
used later on. 



Lemma 5.3. For 1 < k < 2N - 1, there holds 

A fe 



(5.3) 
(5.4) 
(5.5) 



6 



< cosh7 fc 



At 



3 Afc 
2 " T' 



sinh 8k > \ — 



. , OLk / Afc 

sinn — = a / — . 
2 V 10 
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Proof. Invoking (|4.3|) . we have 



cosh 7fc -2 = 7 (^25 + 8X k - 5 



2A* 



V25 + 8A fe 



which immediately implies (|5.3|) . 

The estimate (|5.5|) follows from (|4.1[) by definition. 
Using (|4.3[) . we have 



(5.6) 



cosh7 fc - cosh 4 = 



which together with the definition leads to 

sinh 2 S k = cosh 2 S k - 1 = (cosli7 fc - 1/2) 2 - 1 
= (cosh 7fc + 1/2) (cosh 7fc - 3/2) 
= (cosh 7fc - 3/2 + 2)(cosh7 fe - 3/2). 



Using (|5.3|) . we have 



sinh 2 4 > 2(cosh 7fe - 3/2) > 



This gives (pT4"]) . 

To proceed further, we need the following estimates. 

Lemma 5.4. For 1 < k < 2N — 1, there holds 

k 



exp(-afe) < exp 



exp(-7fe) < ^ ' , exp(-4) < exp 



(5.7) 



(5.8) 

Proof. We only prove (|5.7j) . Other cases are similar. 
Using (|5.5[) and cosha^ > 1, we have 

exp a k = cosh + sinh a k 

2 kn 
= 1 H = sm . 

y/E 47V 



fc 
"5N 



Using Jordan's inequality 



we have 



sinx > —x for x € [0,7r/2], 

7T 



expafc > 1 



VEn' 



For any < x < 2/y/E, we have 

ln(l + x) > x(l - x/2) > x(l - 1/V5) > x/2. 
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Using the fact that k/(v5N) < since 1 < fc < 2N — 1, and combining the 

above two inequalities, we obtain 

exp(-a fe ) < (l + = exp[- ln(l + k/(y/5N))] 



< exp 



2V5N 



sinh a k 



sinh[(M - l)a k ] 
1 



M-l 



□ 



The next lemma concerns the estimate of p k . 
Lemma 5.5. 

(5.9) < 1 - 6p k < I ( ^ + — — - + sinha fc 

6 \ 5 M — 1 

Proof. Using the definition of p kl we get 

l _ 6 = 5 (sinh[A/a fe ] - sinhp/ - l)a k }) 
Pk ~ 5sinh[Ma fe ] + sinh[(M - l)a k ] 

which implies the left hand side of (|5.9I) . Moreover 

5 / sinh[Ma fc ] 

Pfe ~ 6 Vsinh[(M- l)a k ] 
5 

= - (coshafc — 1 + cot[(M — l)afe] sinner) 
Using cosh £ < 1 + sinh £ for any £ G R, we have 
cot[(M — l)a k ] sinhafc < sinha^ 

< sinhafc 

where we have used the elementary inequality 

sinh[M£l 
— A — - > M. 
smh £ 

Combining the above three inequalities, we obtain the right hand side of (|5.9p . □ 



By the definition of p k and the left hand side of ()5.9[) . we get 

(5.10) < Pk < 1/6. 

A direct calculation gives 

(cosh 7 + coshS)((-l) m F m f n - (-l) n F n f m ) 

= Asinh[M7] sinh[M<5] +Scosh[M7] sinh[MJ] 
+ Csinh[Af7] cosh[M<5] + £>cosh[A/7] cosh[M6] 
— sinh[(m — n)8] sinh 7 + (_i) m +™ sinh[(m — 71)7] sinh<5, 
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where 



and 



and 



A = (-l) m cosh[m 7 ] cosh[(n - 1)5} + (-l) m cosh[n<5] cosh[(m - 1) 7 ] 
- (-l) n cosh[n7] cosh[(m- 1)5] - (-1)" cosh[mS] cosh[(re- 1)7], 

B = -{-l) m sinh[m7] cosh[(n - 1)5} - (-1)™ sinh[(m - 1)7] cosh[n5} 
+ (-l)"sinh[n7] cosh[(ra - 1)6] + (-l)"sinh[(n - 1) 7 ] cosh[m<5], 

C = -(-l) m cosh[TO7] sinh[(n - 1)5] - cosh[(m - 1)7] sinh[nd] 

+ (-1)™ cosh[n7] sinh[(m - 1)5] + (-1)" cosh[(n - 1)7] sinh[m<5], 



and 



D = (-l) m sinh[m 7 ] sinh[(n - 1)5} + (-l) m sinh[(m - 1)7] sinh[n<5] 
- (-1)" sinh[n7] sinh[(m — 1)5} — (— l)' 1 sinh[(n — 1)7] sinh[m<5]. 

The following lemma gives a lower bound for \Pkfk — RkPk\- 

Lemma 5.6. There holds 



(5.11) (cosh 7/£ +cosh4)|P fc r fe -i? fe p A .| > — ( 1 - exp 



2M 
'~5N 



exp [M( 7fc + 4)]. 



Proof. A direct calculation gives 

(cosh7 fc + cosh 5 k )(Pkrk ~ RkPk) 
= A k smh[Mjk] sinh[A/4] + B k cosh[Af7 fc ] sinh[M5 fe ] 
+ C k sinh[M 7fc ] cosh[M5 k ] + D k cosh[M 7fc ] cosh[M4] 
+ 2(12 + 2\ k - 72p k ) sinh7 fe sinh<5 fc , 



where 



Ai 



313 



X k + 315 I Pk - 6Xt - 39A fc 



105 
~2~' 



225 

B k = [ 18X k + — + ( ^ + 2X k ) V25 + 8A fc ) Pk sinh 7fc 



229 



25 

— + 4A fc - (14 + 4A fc )V25 + 8A fc ) sinh 7fc , 



C k = -18X k - 



~- 4X k - (14 + 4A fc ) v / 25 + 8A fe J sinh4, 
D k = ((169 + 8X k )pk - (74 + 20A fc )) sinh 7fc sinh<5 fc . 



225 



229 



2A fc V25 + 8A fc p k sinh«J. 
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Using (|5.10|) . we may show 

(5.12) A k ,B k ,C k ,D k < 0. 
Using D k < and invoking (|5 . 10[) once again, we have 

D k cosh[M 7fc ] cosh[M4] + 2(12 + 2A fe - 72p fe ) sinh 7A: sinh 4 

< D k + 2(12 + 2A fc - 72p fe )sinh 7fe sinh 4 

= ((169 + 8X k )p k - (74 + 20A fc ) + 2(12 + 2A fe - 72p fc )) sinh 7fe sinh 4 
= ((24 + 8X k ) Pk - (50 + 16A fc )) sinh 7fc sinh 4 

< - (46 + 44A fc / 3) sinh 7fc sinh S k 

(5.13) <0, 

which together with (|5.12p implies 

(cosh 7fc + coshS k )(P k r k - R k p k ) < 0. 
Combining the above equation with (|5.12[) and (|5.13p . we obtain 

| (cosh 7fc + cosh<5 fc )(P fc r fe - R k p k ) \ = -(cosh 7fc + coshS k )(P k r k - R k p k ) 

= -B k cosh[M 7fe ] sinh[M4] 

- A k sinh[Af 7fc ] smh[M6 k ] - C k sinh[M 7fe ] cosh[M5 fc ] 

+ {-D k cosh[M 7fe ] cosh[M6 k ] - 2(12 + 2A fc - 72p fe ) sinh 7fe sinh 4) 

(5.14) > \B k \ cosh[M 7fc ] sinh[A/4]. 

Using (|5.9[) . we bound \B k \ as 
\B k \ = —B k 

225 /229 A / A 

18A fe + — +(— + 2A fc J V25 + 8A fe J p k sinh 7fe 



25 
T 



4A fc - (14 + 4A fc )v / 25 + 8A fc ^) sinh 



" ( (?2T + T Afe ) ~ X " 7A fc ) sinh 



Ik 

' Ik- 



A direct calculation gives sinh 7 fc > v / 5/2, which together with the above inequality 
yields 

(5.15) \B k \ > | 



By 

exp[2M4] > exp[2kM/(5N)] > exp[2M/(5iV)] 
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It follows from the above inequality that 

smh[MS k ] = i cxp[M4] (1 - exp[-2M4]) 

> i (i _ exp[-2M/(5N)])exp[MS k }. 

Substituting the above inequality and f|5 . 15[) into (|5.14[) implies (|5.11[) . □ 

Next two lemmas concern the upper bounds of Q k and Q m .k- Instead of calculat- 
ing all the coefficients of Q k and Q m , k as we have done for (cosh7fc+cosh5fc)|P^rfc — 
RkPk\, we consider the coefficients of the leading order terms of Q k and Q m ,k- We 
write 

(cosh 7 + cosh<5)((-l) m F m /„ - (-l)"F„/ ro ) 

= -{A + B + C + D)e M ^ +5) 
4 

(5.16) + -{-A -B + C + D)e M ^~ s) + L.O.T. 

with 

A + B + C + D = (e 7 + e 5 ) (f-l) m e -^ n5+m ^ - (_l)™ e -(™7+™5)^ ( 
-A - B + C + D = (e 7 + e -4 ) (_(_i)™ e «<5-™7 + (_ 1 )» e ™-5-«7^ > 

and L.O.T. stands for the terms that are of lower order than e Miri+s ^ and e M ^^ s \ 
Using (|4.14[) . we write (cosh 7fc + cosh 8 k )Q k as 

(cosh 7fc + cosh4)Qfc - Ql cxp[Af( 7fc + 4)] + Q\ exp[M( 7fe - S k )} + L.O.T. 
with 

Ql = i(e 7 « + e Sk ){l2(e? k + e Sk ) + 5{e 2Sk - e 27fc ) - 5(e~ 7fc + e" 5 ") 

— e 1k+Sk (e 7fc + e Sk ) + e 1k ~ Sk — e Sk — 1k \ 
Ql = i(e 7t + e- Sk ){-l2{e lk + e~ Sk ) + 5(e 27fc - e" 2Sk ) + 5(e _7fc + e Sk ) 

We have the following estimate for Q°. 
Lemma 5.7. There exists C such that 
(5-18) \Q k \ + \Ql\<C^X~ k . 



(5.17) 



Proof. We only estimate Ql, and Q\ can be bounded similarly. We firstly write 
Ql as 

Ql = i(e 7fc + e Sk )\[\2 + e- Sk - e 25k )e~< k - (5 + e &k ){e 2 ~< k + e^ k ) 
(5.19) + \2e Sk - 5e~ Sk +5e 25k Y 
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By definition, we have 

(5.20) cosh[27/ c ] = 2 cosh 2 7fc — 1 = cosh^ + \ k + 2, 
which implies 

e 2 7fc _|_ e /fc _ cosh[27fc] + sinh[27fc] + cosh7fc — sinh7fc 
= 2 cosh 7fc + Xk + 2 + sinh j k (2 cosh 7* — 1) 
= 2{e lk + 1) + A fc + sinh7 fc (2cosh 7fe - 3). 

Substituting the above equation into (|5 . 19[) produces 

Ql = i(e 7fc + e'*){(5 - e~< k )[2{e Sk - 1) - e -4 * + e 25fc )] 

- (5 + e 5 * ) [A fc + sinh lk (2 cosh lk - 3)] } . 

Using (|5.6p . we get 

2(e 5,! - l)-e~ 5,t + e 2,5fc = (e 5k - l)(3 + 2cosh4) 
= 2(cosh7fe + l)(e Sk - 1). 
Combining the above two equations, we obtain 

Ql = ( e7fc + e 5fc ){2(5 - e 7fc ) (cosh 7fc + l)(e 5 " - 1) 

(5.21) - (5 + e 5fc )[A fc +sinh7 fe (2cosh7 /! -3)]}. 
Proceeding along the same way that leads to (|5.2ip . we obtain 

Ql = \ (e 7t + e 5k ) { (5 + e- Sk ) [A fe + 2 sinh lk (2 cosh j k - 3)] 

(5.22) + (5 - e»)(l + Se" 5 * + e -M *)(e'* - 1)}. 
Using ()5.6j) gives 

3 

e Afc — 1 = cosh<5fc + sinh 5ft — 1 = cosh7fc — — + sinh^fc. 

Substituting the above equation into (|5.21[) and (|5.22p . and using the esti- 
mates (|53|) and J53]), we obtain ([535]) . □ 

Next we write (cosh7fc + cosh5fe)Q mi fe as 

(cosh7 fc + coshS k )Qm,k = Qm,k exp[M(7ft + 4)] 

+ Q^exp[M( 7fe -4)]+L.0.T. 

with 

Q° m .k = \(e 7k +e 5k ) ((-l) m e- m7fc Qi - e~ mS »Q 2 ) , 
Qi = 3(1 + e 5 ") + e 2S " - e~ 6k - 6p k (5 + e 6k ), 
Q 2 = e 27fc + e-~< k - 3(e 7fc - 1) + Gp k {e~< k - 5), 
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and 

[ Ql.k = + e~ Sk ) H-ire-^Qs - e mS *Q 4 ) , 

Q 3 = 3(1 - 2p)(e- s » - 1) + e- 2S " - e 5 * + 6(1 - 6p k ), 
Qi = -Ca- 

We have the following estimate for Q° ra fc and Q) n k . 
Lemma 5.8. There holds 

1 



(5.23) |Qi| + |Q 2 | + |Q 3 | + |Q4| <C VAfc 



M - 1 



Proof. We only estimate Qi and Q 2 - The terms Q 3 and Q4 can be bounded 
similarly. 

Similar to (|5~20l) . we have 

cosh[24] = — cosh 4 + + 2. 
Using the above equation, we write Qi as 

Qi = 3(1 + cosh 4 + sinh4) - cosh 4 + \ k + 2 + sinh[24] 
- (cosh 4 - sinhtffc) - 6p fc (5 + e &k ) 
= A fc + 5 + cosh 4 + 4sinh4 + sinh[24] - 6p fc (5 + e S ") 
= A fc + (5 + e Sk )(l - 6p fc ) + sinh4(3 + 2 cosh 4). 

Using (|5.20p and proceeding along the same line that leads to the above expres- 
sion of Qi, we obtain 

Q2 = 2 cosh 2 7k + 2 + sinh[27fe] — 3(cosh7fc + sinh7^.) + cosh 7/,. — sinh7fc 
+ 6 Pk (e' yk -5) 
= A fc + 5 - cosh7fc + sinh[27fc] - 4sinli7 fc + 6pfc(e 7fc - 5) 
= Afc + (1 - 6 P fc)(5 - e«>) + 2sinh 7 fc(cosh 7 fc - 3/2). 

Using (gUD and we get 

1 



|Ci| + |Ca| <C V / Afe + 



M — 1 



□ 



To prove Theorem 15.11 we need the following identity that can be found in [5J 
p. 38, formula 1.353(1)]. 

Lemma 5.9. For any g £ (0, 1), we have 

X- ^ n 1 - (g± g3 ) sillx ~ g2jV+1 sin K 2iV + 1^3 + g2jV+3 sin K 2iV ~ 1 

sm^ Ijxj- 1 _9.2 Poe r 9T .i T ,4 



fe=i 



2<> 
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Based on the above estimates, we are ready to prove Theorem 15. II 



Proof of Theorem \5.1\ Using (|5.11|) and (|5.18j) . we have, for m < — 2, 

27V-1 

(5.24) 



\Dy(m,n)\ < ^- ^ exp(-|m|a fe ) sin k 



fe=i 

k odd 



2N 



By (|5~TTj) and ([5T2"5], , we have, for m > 0, 



2N-1 



(5.25) |£>iy(m,n)| < — ^ ( exp(-m7 fe ) + exp(-m<5 fc ) sin — 



fe=i 

A: odd 



(5.26) |D 2 y(m,n)| < — 2J exp(-m 7fc ) sin — + exp(-mfc) sin — 

Z 1 V / 



k=l 
k odd 



Let g = exp[— \m\/(2y/5N)] and a; = 7r/[2iV], and using Lemma 1531 we have 



2JV-1 



k'TT 



N 



J2 exp(-|mK)sin— < J2 ^ 1 sin[(2k - l)x] 



k=l 
k odd 



k=l 



q(1 + q 2 )(1 + q™) 
(1 — g 2 ) 2 + 4g 2 sin 2 x 
2g(l + g) 2 . 

(1 - e 2 ) 2 

2g 



■ SIM 



(1 - 8)'- 



Using Lozarevic's inequality 



coshi < 



/ sinht 



■ SIM. 



and the elementary inequality cosht > e*/2, t <G R, we obtain 

7T 



2p 



(1 - 



■ sin x 
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2 sinh' 



/ |m| 



< 



2iV 



f H \ cosh 2/3 /_H 



2 2 / 3 A0ttN 
< ~3. ex P 



< 



80ttN 



■ exp 
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which together with (|5.24|) leads to(|5.ip. 

For to > 1 and let g = exp[— m/(5N)], we immediately have (|5.2j) . 

The proof of the case when to = — 1 can be done in the same way that leads 



Proceeding along the same line that leads to Theorem 15. 1[ we have the following 
estimate on the solution. 

Corollary 5.10. There exists C such that 



The above estimate suggests that the ghost force actually induces a negligible 
error on the solution, which is as small as e. 
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to (|5.1j) . We leave it to the interested readers. 



□ 



\y(m,n)\ < C- — -, m < -1, 
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